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ABSTRACT 


Nonlinear analysis methods are developed which will enable 
the reliable prediction of the dynamic behavior of the space 
shuttle main engine (SSME) turbopumps in the presence of bear- 
ing clearances and other local nonlinearities. 

A computationally efficient convolution method, based on 
discretized Duhamel and transition matrix integral formula- 
tions, is developed for the transient analysis. In the formu- 
lation, the coupling forces due to the nonlinearities are 
treated as external forces acting on the coupled subsystems. 
Iteration is utilized to determine their magnitudes at each 
time increment. The method is applied to a nonlinear generic 
model of the high pressure oxygen turbopump (HPOTP). As 
compared to the fourth order Runge Kutta numerical integration 
methods, the convolution approach proved to be more accurate 
and highly more efficient. 

For determining the nonlinear, steady-state periodic 
responses, an incremental harmonic balance ( IHB) method was 
also developed. The method was successfully used to determine 
dominantly harmonic and subharmonic ( subsynchronous ) responses 
of the HPOTP generic model with bearing clearances. A reduc- 
tion method similar to the impedance formulation utilized with 
linear systems is used to reduce the housing-rotor models to 
their coordinates at the bearing clearances. 

Recommendations are included for further development of 
the method, for extending the analysis to aperiodic and chaotic 
regimes and for conducting critical parametric studies of the 
nonlinear response of the current SSME turbopumps. 
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NOMENCLATURE 
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Cosine coefficients associated with Master degrees 
of freedom 

Incremental vector of {A^j} 

nth cosine coefficient of the steady state solution 
in x, y-direction at ith node 

Cosine coefficients associated with Slave degrees of 
freedom 

nth since coefficient of the steady state solution 
in x, y-direction at ith node 
Structural damping matrix 

Constant cosine coefficient vector of nonlinear 
force 

nth cosine coefficient of nonlinear force in x, 
y-direction 

Cosine coefficients of nonlinear force associated 
Master degree of freedom 

Cosine coefficients of nonlinear force associated 
Slave degree of freedom 

Generalized damping matrix of rotor model 
Constant sine coefficient vector of nonlinear force 
nth sine coefficient of nonlinear force in x, 
y-direction 

Mass eccentricity of disks 

Nonlinear restoring force vector in x, y-direction 
at ith node 
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= Imbalance forces on rotor 

= The coupling force vector on the housing due to the 
rotor 

= Imbalance force vector of cosine, sine terms 
= Gyroscopic matrix 

= x,y rectilinear housing displacement at ith node 
= Identity matrix 
= Structural stiffness matrix 
= Bearing stiffness 

= Equivalent constant coefficient stiffness matrix 
= Mass matrix 

= Generalized force (modal force) 

= Generalized displacement (modal displacement) 

= Transition matrix 

= Trigonometric coefficients associated with Master 
degrees of freedom 

= Trigonometric coefficients with nth harmonic mode in 
Master degree of freedom 
= Incremental vector form of {Cmn} 

= Trigonometric coefficients associated Slave degrees 
of freedom 

= Overall incremental coefficient vector 
= x,y rectilinear displacement at ith node 
= Total trigonometric coefficient matrix 
= Equivalent trigonometric coefficient matrix 
= Overall coefficient matrix 
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= Time 

= Time increment 
= Time at iT 

r q i 

= {.}/ Generalized coordinates of rotor 

q 

= equivalent trigonometric force vector 
= Overall incremental coefficient vector of nonlinear 
force 

= Constant trigonometric coefficient vector of 
nonlinear force 

= Trigonometric coefficients of nonlinear force 
associated with Master degree of freedom 
= Trigonometric vector of nonlinear force with nth 
harmonic mode in Master degree of freedom 
= Incremental vector of {W^ n } 

= Overall force vector 

= Rational displacement in x ,y-direction at ith node 
= Gap size 
= Subharmonic ratio 
= The rotor state matrix 
= Diagonal eigenvalue matrix 
= Spinning speed of rotor 
= Normalized modal matrix 
= Damping ratio 

= Natural frequency and damped natural frequency 
= Null matrix 
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Superscript 

i = ith node 

y = The node number at which disk is located 

Subscript 

B, b = Bearing related term 
H = Housing related term 
E = External, imbalance 
I = Coupling 
i = ith time increment 

L = Total number of degrees of freedom associated with 
linear motion 

N = Total number of degrees of freedom associated with 
nonlinear coupling nodes 

M = Total number of degrees of freedom (M = L+N) 

R = Rotor related term 
m = Master degree of freedom 
n = nth harmonic mode 
s = Slave degree of freedom 

x = x-direction 

y = y-direction 

x,y,z = x , y , z-direction 
Other Symbols 
{ } = Vector 

[ ] = Square matrix 

P v ] = Diagonal matrix 

{ } = Column matrix 
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[ = Inverse of a matrix 

T 

[ ] = Transpose of a matrix 
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1. INTRODUCTION 


1 . 1 Background 

Modern complex rotating machinery contain various sources 
of strong nonlinearities. These include clearances in bear- 
ings, gears and splines, rubbing in splines, seals and rotor 
blades, squeeze film dampers and other fluid effects. Observed 
nonlinear behavior of actual rotor systems includes jump 
discontinuities (Erich, 1966), large subsynchronous motion 
(Bently, 1979, Muszynska, 1984, and Beatty, 1985), quasi-peri- 
odic response and possible chaos (Neilson and Barr, 1988). As 
stated by Nataraj and Nelson (1987), the future developments in 
modern machines heavily depends on the ability to identify, 
understand, mathematically model and analyze systems involving 
nonlinear components. 

The turbopumps of the space shuttle main engine (SSME) are 
complex rotor-housing systems. Their response and stability 
can be critically dependent on the coupling forces between 
their flexible rotors and housings, transmitted through the 
working fluids, seals and bearings. Of particular concern is 
the effect of the existing essential deadband clearances 
between the ball bearings and housing. These clearances cause 
a turbopump rotor to respond as a nonlinear dynamic system. In 
general, the rotor-housing system could, consequently, exhibit 
undesirable subharmonic, combination and internal resonances 
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(Tondl, 1965). It could also respond in an aperiodic or 
unpredictable chaotic fashion (Shaw and Holmes, 1983 and 
Neilson and Barr, 1988). The transient response of the rotor 
in passing through a critical speed can be quite different in 
the presence of a nonlinearity ( Ishida et. al . , 1987). 

Childs and Moyer (1984), using a Runge-Kutta integration 
scheme, carried out a transient nonlinear analysis of the HPOTP 
(high pressure oxygen turbopump) of the SSME and reproduced a 
subsynchronous rotor response similar to that observed in 
experimental test results. They attributed the subsynchronous 
response to the presence of the bearing clearances and desta- 
bilizing impeller-diffuser forces. Glease and Bukley (1984) 
examined the effects of bearing deadbands on a modified Jeff- 
cott model subjected to seal cross coupling forces and 
described responses of the rotor to imbalance forces of the 
subsynchronous, synchronous and nonperiodic type. Using a 
similar model as that of Glease (1984), Day (1987) and Zalik 
(1987) addressed some aspects of the response in the presence 
of bearing clearances and cross coupled seal forces as leading 
to a limit cycle with a "nonlinear natural frequency." For 
certain ranges of the systems parameters, quasi-periodic 
response of the rotor can occur. 

Quite often, it is essential to determine the steady state 
periodic response of rotor systems in the form of self excited 
limit cycles or forced motion due to rotating unbalance. 
Accurate prediction of the nonlinear periodic responses and 
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their stability plays a central role in developing a complete 
picture of the dynamic behavior of nonlinear rotor systems as a 
function of their parameters. Efficient general methods are 
obviously needed for determining directly possible periodic 
solutions of large rotor-housing systems involving nonlineari- 
ties. In contrast, a direct numerical integration method used 
to obtain the steady state periodic solutions: 

a. Will require much longer computational time for the rela- 
tively lightly damped rotor systems to settle to a steady 
state . 

b. Will determine only stable periodic responses and not the 
unstable solutions which may be used to indicate the 
onsets of transition to chaotic responses. 

c. Might miss a close-by large steady state solution from 
among multiple nonlinear solutions due to an unfortunate 
selection of initial conditions. 

1 . 2 Objectives and Outline of Study 

The main objective of this study is to develop reliable 
and efficient analytical-computational methods of the nonlinear 
dynamic analysis of large, rotor-housing systems such as the 
turbopumps of the space shuttle main engines (SSME). 

The present study examines some aspects of the nonlinear 
behavior of a general multi-disk rotor-housing system. For 
that purpose, a generic model of the SSME turbopumps proposed 
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by Bavis et al. (1984) is used to which new nonlinear analysis 
methods developed by the present author and co-researchers are 
applied. One of the methods extends an incremental harmonic 
balance method developed earlier by Choi and Noah (1987) to 
the determination of the periodic response of multi-degree of 
freedom rotor systems. Associated with this method, a Quasi- 
Newton method was also developed to enable the determination of 
possible multiple solutions for a given set of rotor parame- 
ters. For the transient analysis, a convolution method is 
fully developed which takes advantage of the fact that a turbo- 
pump system consists of a linear rotor model coupled to a 
linear housing model through a local nonlinearity, that is, 
bearings with clearances, to construct a highly efficient and 
accurate computational procedure. The method can be equally 
applied in presence of other coupling nonlinearities between 
the rotor and housing. 
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2. NONLINEAR ROTCRD YNAM I C ANALYSIS OF LARGE SYSTEMS 
2 . 1 Transient Analysis Methods 

The direct numerical integration for determining the 
transient response of larger nonlinear rotor-flexible housing 
systems may require excessive computational time and involve 
unacceptable round-off errors. To remedy these problems 
different procedures have been developed by analysts to deter- 
mine the transient response of large order rotor systems. The 
procedures can be recognized as falling under one of two basic 
approaches. Those using physical or modal coordinates of the 
complete system and those using the coordinates of the individ- 
ual components of the system. The methods also differ in the 
numerical integration methods selected for the analysis. A 
comprehensive review of existing transient analysis techniques 
can be found in a report by Noah (1986). 

An effective method which proves to be highly efficient in 
determining the transient response of linear systems under 
specified general excitations is by utilizing the Duhamel 
integral/ or more generally/ the transition matrix of the 
system (Meirovitch/ 1980). Von Pragenau (1981) utilized the 
transition matrix with a linear system, stating that it offers 
the simplicity of the Euler method without requiring small 
steps. Von Pragenau maintained that for systems with constant 
coefficients, the stability and accuracy of the method are 
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acquired through the closed form solution of the transition 
matrix. 

In the present study, the convolution methods ( Duhamel and 
transition) are shown to be very effective when extended for 
application to linear systems with local nonlinearities . The 
forces at the nonlinear locations are treated as the external 
forces on the systems or subsystems and iteration is used at 
each time increment to determine the magnitude of the forces 
for subsequent increments. 

The convolution approach is applied to a general rotor- 
housing system with rotor imbalance during start-up or shut- 
down. In the present work, eigen-coordinates are used to 
represent both housing and rotor. The local nonlinearities 
were taken as deadband clearances at the rolling element bear- 
ings which support the rotor in its housing as in the SSME 
turbopumps. The integral formulation of the rotor motion is 
represented by its transition matrix and that of the housing by 
a Duhamel integral. Kubomura (1985) used a Duhamel integration 
method to achieve dynamic condensation of a substructure to its 
coupling points to other structures. Convolution was also used 
by Tongue and Dowell (1983), and Clough and Wilson (1979) to 
reduce system coordinates to that of the nonlinearities. 

2.2 Periodic Solutions 


Several methods have recently been advanced for deter 


mining the periodic response of low order nonlinear rotor 
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systems (Yamamoto 1954, Childs 1982, Saito 1985, and Choi and 
Noah 1987). For application to larger, multi-disk rotor sys- 
tems, Nataraj and Nelson (1987) developed a periodic solution 
method based on a collocation approach for the response of the 
rotor. They utilized a subsystem approach (Cipra and Uicker, 
1981) to reduce the size of the resulting system of algebraic 
equations . 

In this report, a general approach for the determination 
of the synchronous and subsynchronous response of complex rotor 
systems is presented. The present work extends the harmonic 
balance method developed earlier (Saito 1985, and Choi and 
Noah 1987 and 1988) for a modified Jeffcott rotor model with 
bearing clearances to the determination of the periodic 
response of nonlinear multi-disk rotor systems with flexible 
housings . 

It is shown that the harmonic balance method, as applied 
to a large, rotor-housing system with bearing clearances, can 
be made to be highly efficient. This is achieved by using a 
version of an impedance formulation in which the system is 
reduced to its displacements at the bearing clearances. (See 
Noah, 1984 and Fan and Noah, 1989). In this case, the 
impedance method is applied to each of the harmonic components 
of the assumed periodic solution. 
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3. APPLICATION TO ROTOR-HOUSING SYSTEMS 


3 . 1 The Model 


A representative complex rotor system with flexible 
housing is shown in Fig. 1. The particular model shown in 
which the present method can be readily applied/ represents the 
high pressure oxygen turbopump (HPOTP) of the space shuttle 
main engine (SSME). The interaction forces between rotor and 
housing include various seal, impeller, turbine tip clearance, 
bearing clearance and fluid side forces. 

a. - Rotor 

The equation of small transverse motion, {r} of the rotor 
under transient external and imbalance forces may be written as 
[M] r {R} - i [G] {R} + [K] r {R} = { F x } R + {F E } R (1) 

In equation (1), {r} represents the translational and rota- 
tional displacements and only includes those displacements 
which are associated with masses and diametral moments of 
inertia. Other displacements are reduced out using static 
condensation. The spinning speed of the rotor is denoted by 
(see Fig. 2), [G] is the gyroscopic matrix corresponding to the 

attached disks, { f j} r represents the coupling forces on the 
rotor due to coupling to the housing while { F R } R represents the 
external forces including the imbalance forces. 


Let {R} = m R { q} R 


( 2 ) 
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in which ['Hp is the modal matrix, normalized with respect to 
the mass matrix of the nonspinning rotor and {q} n is the asso- 
ciated modal coordinates. 

Using the modal transformation (2), equation (1) is 
written in the form 


iq) R - 


m R * [gi d m D { q } D + pAu 


R 




({ F I^ R + i F e^R ) 


(3) 


Where P A v] R = 


“ n are the natural frequencies of the nonspinning of the rotor. 
Adding a modal damping matrix PCv] to equation (3) yields 

iq] R - [D] r (q] R ♦ PA.] r lq} R = (*]* - U E J R > <41 

where [D] R = fC.] R - [*]£ i [G ] R [f] R (5) 

PC.] R = P 2Co) n .] R (6) 

b. Housing 

The equation of motion of the housing may be expressed as 
[M] h { H} + [K ] H { H} = {F X } H (7) 

where {h} represents the transverse and rotational displace- 
ments at each node in the Y and Z directions. In terms of the 
modal coordinates of the housing while uncoupled to the rotor, 
equation (7) takes the form (after adding a modal damping 
matrix) 
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{q} H + [di h i q} H + {q} H = } H + {f e } h ) (8) 

Where [D] H = ['CJ H = t'25w n J H (9) 

['A.1 h = ['«n-l H (10) 

c . Coupling Forces 

Forces due to linear coupling between rotor and housing 
may include the direct and cross coupling forces due to seal, 
impeller and turbine forces. The nonlinear coupling forces are 
taken as those due to the clearances between the rolling 
element bearings outer races and the housing. Fig. 3 shows a 
model for the gap at each of the loosely supported bearing. 

The bearing force acting on the housing due to the 
stiffness, shown in Fig. 4 in the Y-direction, is 

(F B ) Y = (R " 6 } — - Hy for R > 6 

R 

( F fi ) y = 0 for R _< <5 

Analogous equations hold for the Z direction, in which K is 
the stiffness at the bearing support and 

R = 7 (R y - H y ) 2 + (R z - H z ) 2 






Figure 


Rotor and housing displacements at bearing location. 
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where 


Ry is the physical displacement of rotor in 

Y direction at the bearing location. 

H y is the physical displacement of housing in 

Y direction at the bearing location. 

Rz is the physical displacement of rotor in 
Z direction at the bearing location. 

Hz is the physical displacement of housing in 
Z direction at the bearing location. 


The bearing forces in the Y and Z-directions can be written as 


( ^ y ^BY ^ ^Y 


( f b j z k bz (r z 


Where 


k b = k b (1 


K 


B 


h y ) 


V 


P f ° r 


for 


R _> 6 
R <6 


(ID 

( 12 ) 

( 13 ) 


0 


( 14 ) 
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3 . 2 Method for Transient Analysis 


A hybrid method proposed by Noah e_t al_. ( 1986) and (1988) 
is utilized here for the analysis of the coupled rotor-housing 
system. The displacements of the housing are best represented 
by a Duhamel integral. Due to the presence of the skew-sym- 
metric gyroscopic terms, the equations of motion for the rotor 
are transformed to first order. The displacements of the rotor 
are expressed in terms of the transition matrix for the motion 
of the rotor. 


3 . 2a Transition matrix for the rotor 


The equations of motion of the rotor, equation (4), are 
cast in first order form, 


1 ^ 1 j> = f r { u 1 r + Ip. 


where [ a] p = 


[0] ['IJ 

- [D]„ 


(16) 



0 

t R ( { F I + { F E 


(17) 


If the spinning speed is a function of time, the gyro- 
scopic term of equation (5) is moved to the R.H.S. of equation 
(4). In this case, the coefficient matrices, given by equa- 
tions (16) and (17), are replaced by 
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to] no 

■ pAs] R - rC ' ] R 


( 16 ' ) 


and 


0 


0 

^R ( ^ F I + ^ F E ^R ^ 

+ 

1 1 

m* hG] E m R tq} R 


(17' ) 


In that case it may be computationally more efficient to use 
the Duhamel integral, rather than the transition matrix formu- 
lations to represent the rotor. The formulation using the 
Duhamel integral is applied to the housing as is shown later. 

In general, however, if the physical rather than the modal 
coordinates are used to describe the motion of the rotor, or if 
$ is constant and is kept at the left hand side of the equa- 
tions of motion, the transition matrix and not the Duhamel 
formulation method must be used. 

The solution of equation (15) can be written as 


[«]»t. , .t E«] «( t-T ) 

R = 6 ‘'O' 


lu(t)} T , = e* JR {uj + / e 1 Jp '{P(t)} p dT 


( 18) 


where 


[a] R t 


k=0 


[«]* t k 
k! 


= transition matrix (19) 


and { u 0 } is the initial generalized coordinates 


Equation (18) can be cast in discretized form in which the 
generalized forces, { P } _ , are taken as varying linearly between 
time ti and ti + i so that 
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fc - fc i 

{P(t) } = {P(ti)} + ■ - - T - - ■ ({P(t i+ i)} - {P(ti)}), 

ti < t < t i+ i 
and ti + l = ti + T 

Where T is a small increment in time. The discretized 
form of equation (18) can then be shown to take the form 

t U(t i + l>} R - l u<t i>} R + - l' 1 -" Ul^dPIti)) 


+ M 


-1 {^WIr - { p(t i ) }r , r ,-1 


R 


) - [ ° ] R U P(t i + lU R - { p(t i> }r) 

( 20 ) 


where [Q(T)J = e 


‘«Jr t 


( 21 ) 


3 . 2b Duhamel integral for housing 

The symmetric form of the equations of motion of the hous- 
ing allows expressing the modal displacements of the housing by 
a Duhamel integral. The integral formulation has the advantage 
of providing a closed form expression as opposed to conven- 
tional numerical integration schemes. Also this representation 
allows dealing directly with diagonal matrices which results in 
higher computational speed and accuracy. 

Based on equation (8), the housing generalized coordinates 


can be expressed as 
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-£co n t 

{q(t)} H = ['e cos w d tj {q(0)} H 

, — • 

+ r-^i e sin u> d tj ({q(0)} H + { < 3 ( 0 ) } R ) 

t , -5co n (t-t) 

+ / ['— e sin ui d (t-x) J { P( x ) } dx (22) 

o d 

• 

where { q ( 0 ) } , { q ( 0 ) } are initial solutions, u> n , <u d are 

undamped and damped natural frequencies of housing, 

£ is damping ratio of housing, { P ( t ) } H = ['F]g({F I } H + 
{F } H ) the generalized forces act on housing. 


The expression (22) is next written in discretized form, 
with the generalized forces taken as before as varying linearly 
within each time step, or 

-ati+l 

{ q ( ti + i ) } H = ['e cos bti + iJ {q(0)} H 


1 ^^ 1+1 • 

+ P^ e sin bti+i^ ] ( {q( 0) } H + a{c(0)} H ) 


sin bt . . cos bt. . 

; 2 ' 2 * EA i + l + 2 — • EB i + l } 

b(a Z +b Z ) 11 b(a z +b z ) 


( 23 ) 
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{q(t i+ i>} 


H 


r-e 


“ at i+l 


( a cos 


bti + i + b sin bti + i). ] {q(0) } H 


+ [ v e ati + 1 (- sin bti + i + cos bti + i ) . ] ( { q ( 0 ) } R + a{q(0)} H ] 


H- 


-a sin bt i+1 + b cos bt i+1 _ 

+ { 22 * EA i + l 

b(a^+b^ ) 


a cos bt i+1 + b sin bt. +1 } 

+ s 5 . EBi + 1 

b(a z +b z ) 


(24) 


a 

= 


(25) 

b 

= 

<*>d 

(26) 

EAi + i 

= 

e a ^. EAi + A i+l 

(27) 

EBi + i 

= 

© . EBi ®i + l 

(28) 

EAi 

= 

e “ at i-l Al + e " ati “ 2 A 2 + + Ai 

(29) 

EBi 

= 

e -at i-l Bi + e" ati ~ 2 B 2 + + Bi 

(30) 

A i + 1 

= 

Pi +1 (a cos bti+i + b sin bti+i) 


p . 


p . 



i- (a 2 cos btf + i + 2ab sin bti + i - b 2 cos bti+i) 

T(a 2 +b 2 ) 


+ — — - e~ aT (a 2 cos btf + 2ab sin bti - b 2 cos bti) 

T(a 2 +b 2 ) 


(a cos bti + b sin bti) 


(31) 
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B-[ + i = Pj. + i (a sin bti + i - b cos bt^ + i) 

P i+1 " P i 2 2 

- x — x — (a sin bti+i - 2ab cos bt^ + i- b sin bt^ + i) 

T(a z +b^) 

P . - p . 

+ 1 + ~ — 5— ^ e” aT (a 2 sin bt^ - 2ab cos bti~ b 2 sin bt^) 

T(a +b ) 


-aT 


- Pi e (a sin bt^ - b cos bt^) 
P is element of vector {p}. 

3 . 2 c The coupling forces 


( 32 ) 


The coupling forces acting on the housing are given by 

l F I } H = CK] x ({R} - {H}) + [C] T ( {R } - 1 H } ) 
and those on the rotor are 


( 33 ) 


! f i1r " - < F iJ 


H 


Where [K]j and [CJi are the coupling stiffness and damping 
matrices, respectively. The coupling stiffness matrix includes 
the updated bilinear bearing stiffness at each iteration. 

3 . 2 d The computational procedure 

For the responses at time t = the coupling forces 

between the housing and rotor are unknown. An iterative tech- 
nique is used to calculate these forces. The following itera- 
tive loop is used in both the Hybrid and Runge Kutta methods: 
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If the responses at t = are known and responses at 
t = ti+i are desired then: 

(a) Set { Fj ( t i+x ) } = { Fj( ti) } 

(b) Calculate {p}r and { p} jj 

(c) Solve {u(ti+i)} R and {q(ti+l)} H , {q(ti+i)} H by the Hybrid 
method (or the Runge Kutta method). Calculate the 
Euclidean norm of {u(ti+i)}R and {q(ti + l)}R. 

(d) Transfer {u(ti+i)} R and {q(ti+i)} H , {q(ti+i)} H back to 
physical coordinates, then a new connection force 

{ F i ( t i+i ) } can be calculated. 

(e) If the Euclidean norms of two consecutive iterations of 
both rotor and housing are close enough then stop the 
iteration and go to step (f), otherwise go to (b). 

(f) Move a time increment and go to (a) until t reaches a 
specified time for terminating the run. 

3 . 3 Method for Obtaining Periodic Responses 

A typical multi-disk rotor interacting with its flexible 
housing through rolling element bearings with deadband clear- 
ances is considered. The equations of motion of the rotor can 
be recast from equation (1) in matrix form as 

[ M ] p { R} + [ C ] p { r} + (K ] R { P} = { F c Jr + ^ F i^r (34) 

Where [C]r includes damping and gyroscopic terms. The 
displacement vector {r[ is defined as 
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{ R 1 " l R lx' R ly ' * * ’ '&My /6 MxJ 

where and gi stand for displacement and rotation, respec- 
tively, and subscript M represents the total node number. 

In particular, the forces in the y-direction on the bear- 
ings with clearances are expressed as 

F bmy = K bmy ( R my " H ny ) d ~ ■ ) (35) 

y ( R my~ H ny) + (%iz~ H nz) 

where m and n denote the bearing node numbers on the rotor and 
housing, respectively, and y and z denote the y-axis and z-axis 
of the rotor and housing coordinates. The rotor displacements 
are denoted by R while H denotes the housing response. The 
symbol 6 represents the size of the deadband clearance between 
a bearing and the housing. 

The bearing forces, F^y or Fb mz , will vanish if 6 is 
larger than the radial displacement of the rotor relative to 
the housing which is the denominator in equation (35), other- 
wise the bearing forces will be as given by equation (35). A 
similar expression can be written for the nonlinear bearing 
forces in the z-direction, which are denoted by Fb mz . 

The equations of motion for the housing can be written in 
the matrix form as 

[M] h {Sj + tc] H { H } + [K] h (K} = (Fj) h (36) 

where [CJh represents proportional equivalent viscous 
damping coefficients. 
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In order to reduce the system to its displacements at the 
nonlinear bearing supports, the equations of motion need to be 
modified. This is achieved using a subsystem approach. The 
coordinate vector can be modified so that the transverse 
displacements at the gap are listed first 

1 R } = { R ly , 2y , . . . ,Ny ' R lz , 2z , . . . ,Nz ,R ( N+l ) y , (N+2 ) y , . . . ,M y ' 

R (N+1) z ,(N+ 2) 2 ,...,M z } T 

where (1,2,... ,N) is the bearing node numbers related to non- 
linearity while (N+l, N+2, . . . ,M) is the other node numbers 
involving linear displacements which will be eliminated in the 
assembling procedure. 

The housing coordinate vector can be rearranged in a simi- 
lar fashion. After rearranging the coordinate vectors, the 
matrices [M] R , [M] , [C] R , [C] H , [K] R and [K] H , and other force 
vectors also need to be modified according to their vector 
components . 

3 . 3a Computational harmonic balance method 

Extending the procedure developed by Choi and Noah (1987), 
a steady state, a periodic solution for the motion of the rotor 
can be represented by a Fourier series expansion. For the 
displacement of the ith node, one writes 
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R iy ' A 0y + 


I [A* cos ^ sin 


n=l 


ny 


ny 


(37a) 


Riz — + 


Oz 


" ( A n Z cos + B n, •!" 


n=l 


nz 


(37b) 


where v is the subharmonic ratio, which is unity for harmonic 
and superharmonic cases, or an integer for subharmonic cases. 
Similar equations are written for the displacements H i y and Hj_ 2 
of the housing. Since the motion is periodic and steady state, 
the nonlinear coupling forces can also be written as 



N 

z 


n=l 



cos 


ntot 

v 


• 1Z 



N 

E 

n=l 



cos 


nwt 

v 


+ D 


ny 


+ D‘ 


nz 


sin 


ntut \ 
v ' 


sin 


ncut 

v 


(38a) 


( 38b) 


The nonlinear forces exerted on the housing are of equal and 
opposite sign to those of equations (38). To apply the har- 
monic balance method, one forms the vectors of the coefficients 
of cos — — and sin — — , for n=l,2,...,N. 


The constant coefficients of the displacements of the 
coupling forces are then written as follows: 

For the constant term, 


IKJr 




(39) 


where 
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- 

(a 1 a 1 a n a n 1 t 

K°i ■ 

j a n+1 a n+1 a m a m I t 

1A 0y ' Oz ' * * ’ ,A 0y ,A 0z ; 

- 

{ r ^ c N C N 

tL 0y' L 0z' * * * ,L 0y '^Oz ‘ 

lc s °] - 

rn 

{ 0 , 0 , 0 , • • • / 0 f 0 } 


and the subscripts of "m" and "s" stand for the master and 
slave part for a reduced algebraic system resulting from apply- 
ing the harmonic balance method. The reduced out slave coordi- 
nates correspond to those coordinates where no nonlinear coupling 
forces exist. Therefore, the stiffness matrix of the rotor, [K] R , 
is partitioned as 


tK1 F. 


, ( 2Nx2N ) 
mm R 

[K 

[K 1^ 2LX2N) 

[K 


ms' 


SS' 


(2Nx2L) 

R 

( 2Lx2L) 
R 


where N is the total number of degrees of freedom at the 
nonlinearity and L is the total number of degrees of freedom 
associated with linear displacements of the rotor. Equation 
( 39 ) is expanded as 


[K ] {AM + [K 
mm R m 


ms 





(40a) 


[K ] {A^} 

sm R m 






(40b) 
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Combining equations (39) and (40), the following equation is 
obtained 


['We ' [K ms>R lK ss ! ; llK sn, ! pH^l = t C “l 


m 


(41) 


By setting, 


‘ K >P. ■ "W* - f K ms'R [K s S J; 1 [K s m ’ R ' (2H x 2N) 


the final reduced matrix can be written as 




(42) 


Because of the bearing clearances relations already given by 
equation (35), a discrete FFT algorithm will result in {C^} 
being a function of all the Fourier coefficients in equations 
(37) and the corresponding coefficients for the housing 
displacements. The implicit equation form of the relation (42) 
can be then expressed as 






a n a n 

* ' A ny ' nz 


)} 


(43) 


where N is the total nonlinear degree of freedom number and n 
represents the retained harmonic terms. Newton-Rapshon (2nd 
order) iterative method is used, and for that purpose an incre- 
mental form is written in place of equation (42). 

% 

If one sets 


{ A U } = {A U } + { AA } 

1 m J 1 m 1 1 m J 


{C° } = {C° } + {AC } 

1 m J L m J 1 m J 
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then 

<*1r 1 4 V - t iC m) " (O - [K1 Rt A mJ (44) 

This equation is solved sequentially with an equivalent set for 
the housing at each increment and iteration is used to obtain 
compatible coefficient increments for the rotor and housing. 

Applying the harmonic balance method to the equations of 
motion (34), the cosine terms lead to 

- (™) 2 [M] r {a} + (^)[C] R {B} + [K] R { A} = { C n } + { g c } (45) 


where 

{A} 

{ B} 





B' 


nz 


. /A' 


M 

ny 


. t B 


M 

ny 



T 


T 


(M = N+L) 


1 C n J ' l c Jy' C L' ,, * #C ny' C nz' 0 ' 0 ' 


,N _N 


{ g c } — ! 0,0,0 / »..,( mew )"^ " , 0 , ( mew ) ^,0,»..,0} 

if n = v (n=l, 2, . . . ,N) 


= { 0} if n / v 

where the y's are the nodes at which the disks are located 
and M is the total node numbers and N is the node on which the 
nonlinear coupling forces of the rotor to the housing exist. 
Similar expressions are written for the housing. 
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For the sine terms, one obtains 


- (^-) 2 [M] C {B} - <™)[C] r {a} + [ K ] R { B } = {dJ + {g J 


n J 


where 


{ D } 
1 n J 


l D ny' E n 2 ' 


0 n l ^ 

* ' ny ,L nz' u ' * * * ' u ' 


{ g g } = |0,0,0,*»*,( me ) , 0 , ( me (Ai ) 

if n=v (n=l,2, . . . ,N) 

= { 0} if n/v 

Combining equations (45) and (46) 

[S] R { Q} = { w} + {u} 

where 

i T7 1 J /iN /-iN n N _.N n nl T 

t"l = t C ny' C n z' D ny' D nz' 0 '-"' 0 l 


to) = (A N ,a n ,B N ,B K B M ,B H } T 

1 J L ny nz ny nz ny nz J 


To apply the reduction technique, equation (47) can be 
tioned as follows. 


[S]< 4Nx4N) 
mm R 

rc , ( 4Lx4N ) 
sm R 


r o -j ( 4Nx4L ) 
S ms R 

rc , ( 4 Lx 4 L ) 
lS ss J R 


/ ^ °mh 
‘(Osl 5 


,W, , ,t U m) 
1 o J 1 M " S 1 


where 




1 A 1 

,b n ,b n } 

ny nz' 

ny nz J 

N+l ,N+1 

b m 

, . • * , O , 

ny nz 

ny 


(46) 


(47) 


parti- 


) (48) 
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{W } = {c 1 fC 1 , . . . ,d n ,d n } t 

1 m J 1 ny nz ny nz J 

{U s } = {0,0,...,(mea J 2 ) Yl / 0,0,(mea ) 2 ) Yj / 0,...,0} T 

{ Ujjj } = {0,0, •• .,0,0} 

Equation (48) can be expanded as 

‘WV + ‘WV * (V + IV 
'WV + ‘WV ' IV 

From equation (49b) 

«V - ^ss'r 1 t U s 1 - ‘ S ss J R 1[S ms)Rl°ml 

Combining equations (49a) and (50), the following equation 
obtained. 

I'Vr ' s n s J R tS ss ] R llS .m I RlV = IV + IV 

- ‘VW^lVl 

Set 

tS>R * »Wr ' W«,JRW..^ 1 f 8 «.>R' < 4N X 4N) 

(“1 ■ IV - ‘WVV'lV' (4N * 11 

Equation (48) can then be rearranged as 


(49a) 

(49b) 

(50) 
is 

(51) 


is’rIV - tv + ( U 1 


(52) 
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To apply a numerical iteration procedure, one can use the 
following incremental form 

1C } = {Q m > + {*Q 1 
mm m 


tw,,,} = (W m } + tiMj (53) 

Substituting equation (53) into (52) 

[S] R {Ao } - {AW } = {w } + { U } - [S] {q } (54) 

R m mm k m 

Since the equation (54) is obtained only for one harmonic term, 
one can expand equation (54) for general retained harmonic 
terms 


‘V 


R 


{Ac } " t A W } = {w } + {u } - 


mn 


mn 


mn 


n 




le 

R lv mn 




(55) 


where n = 1,2, . . . , N . 

Equation (55) represents the cosine and sine incremental 

terms and equation (44) represents the incremental constant 

terms, so one can combine these equations. {q } and {a } 

mn m 

contain the final solution forms which are solved by the 
Newton-Rapshon method. The overall incremental form for the 
rotor system will be as follows: 

[T] R {A r } - {Av} = { z } (56) 


{Ar } 






A A N a a n 
A 0y ' Oz 


f 


where 
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and 


AA ly' AA lz' AB ly' AB l z ' * * * ,AB ly' AB lz' 


AA ny' AA nz' AB ny' AB nz' * * * ' AB ny' AB nzJ T 


{Av} ~ { AC 0 y' CA 0z' AC 0y' AC 0z' * * *' AC 0y ,AC 0z' 


AC ly ,AC lz ,AD ly ,aD 1z ' * ’ * ,AD ly , aD 1 Z ' 


AC Jy' AC nz' AD iy' AD iz AD nv' AD n*l T 


ny rw nz^ 


tT] R = 


[K]. 

o J 

0 


[Si] 


• • • • 


0 

0 


[ s ] p 

n K 


where 


I 2 ) - ti C °} - '{ W ml I - tUj} - ( S ll P tC ml ). 

i W m2^ i U 2^ " [S 2 ] r{ C m2l ' * * * ' 
l W mn) - 1*3 - lV R («m„l] T 

[K] p is ( 2N x 2N ) and [ST] p is ( 4N x 4N) for (i = l,2,...,n). 
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From equation (56), { Av } is written with partial deriva- 
tives of {^r}. The partial derivative terms yield the Jacobian 
matrix, so one can finally obtain the following equation 


{Av} = [p] { r}’ 


(57) 


where 


[P] = 


3 c 1 

3A 1 

Oy 

3C 1 

3aJ 

Oz 

3C 1 

_I0z 

3A n 

Oy 

3c 1 

3b N 

nz 

3C 1 

C 0z 
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C 0z 
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L 0z 
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Oz 
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•** 3b N 

nz 

• 

• 

3d" 

nz 

• 

3 D N 
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3D* 1 

nz 
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3 C N 

nz 
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3A i 

Oz 
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• • • XT 

»b n 

nz 


where the superscript stands for the degree of freedom number 
while subscript stands for the retained harmonic term in y or z 
direction . 

The [P] matrix can be determined using numerical differen- 
tation of the discrete and inverse discrete FFT method. The 
matrix has to be updated at each iteration step until it 
converges. The steady state solution procedure for flexible 
housing is the same as for the rotor system as outlined above. 

After calculating the rotor response using an iteration 
procedure, the coupling between rotor and housing are auto- 
matically calculated. These nonlinear coupling forces are then 
applied to the housing system to update the absolute housing 
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displacements and nonlinear coupling forces. Another iteration 
process is used to obtain once more the updated rotor response 
and then iterate until both rotor and housing responses have 
the same steady state convergent values. 

3 . 3b Determination of multiple solutions; A global 
Newton's method 

A method is further developed and adapted to rotor systems 
which is capable of locating all possible multiple solutions 
for the rotor's response. The method uses a double dogleg 
scheme (Dennis and Schnabel, 1979) with a local minimum 
algorithm to determine the solutions. A tunnelling method, 
developed by Levy et al. (1978) is further developed and used 
to remove the solutions already found from the formulation. 
Preliminary testing of the method showed the method to be 
capable of obtaining multiple solutions of the Duffing type 
equations . 

An exposition of this method and preliminary results on 
its application will be published in due time. 
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4. NUMERICAL EXAMPLES AND DISCUSSION 
4 . 1 Transient Response 

The present convolution method is applied to a modified 
version of a rotor model constructed by Davis et al. (1984). 

The model was proposed to represent a simplified generic model 
of the space shuttle main engine turbopumps. The parameters 
and coefficients of the present generic model (shown in Fig. 5) 
are given in Tables 1 and 2. The imbalance forces are taken as 
shown in Fig. 5. 

Response of the generic rotor-housing model at the bearing 
location is determined for the hypothetical start-up-shutdown 
case shown in Fig. 6. For the start-up case, from 0 to. 0.01 
seconds, using an increment of 1x10"^ sec., a comparison was 
made of the computer CPU time (on a VAX 8300) for the Runge-Kutta 
4th order and the present convolution method (Fig. 7) . It is 
seen that the hybrid convolution method is approximately twice as 
fast as the 4th order Runge-Kutta method. In addition, a desired 
accuracy can be met' by the convolution method using a much 
larger time step than that required by the Runge-Kutta method. 

Figs. 8-a and 8-b show the radial bearing force in bearing 
1 for the linear and nonlinear (in presence of a clearance) 
cases, respectively, for the fast start-up of Fig. 6, 0-0.3 
seconds. The initial conditions were taken as zero. 

For the nonlinear cases, the radial force in bearing 1 
has a peak close to <j> = 2600 rad/sec. (24830 rpm) during 
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ROTOR: Shaft diameter 00 ■ 3.0 inches, 10 * 0.0 inches 
Material: Steel, E - 3.0 x 1 0® psi 
Joint length ■ 3.0 inches 
Rotor length > 27.0 inches 

HOUSING: Housing/Rotor weight ratio ■ 6/1 

(EI) r - 1.1928 xIO 8 lb-in 3 
< e,) h - 2.4 x10 s lb-in 3 

F^tJ-m^e^cos^t Kg^* 4.0x10* ih/in 

F 2 <t) ■ m^e ^ sin p t Kgy^- 5,8 * ^ ,n 

Fj(t) m mj® p 2 cos $ t 
F 4 (t) • n^ p 2 sin p t 
Fg{t) m mge p 2 cos 1 1 Kp 7 ■ Kgy m 5.0 x 1 0® Ih/in 
F 6 (t)-m 3 ei 2 sin*t 6 = 5.0x10“* in 


Figure 5. The generic model. 


Kg^-iOxItr ItVin 
Kgy - 1.5x10 s Ih/in 
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Table 1. Parameters of the Generic Model 



Maas 

Imbalance 

Moment of Inertia (lb-ln-sec 2 ) 


(lb-sac 2 /In) 

(In) 

Polar 

Diametrical 

Disk 1 

0.0259 

0.002 

0.028 

0.0159 

Disk 2 

0.0389 

0.002 

0.193 

0.101 

Disk 3 

0.0518 

0.002 

1.0057 

0.5078 


Table 2. Coefficients of the Generic Model 
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start-up and a much higher peak at 4> = 3560 rad/sec. ( 34000 
rpm) during the (slower) shutdown in Fig. 9. The linear system 
with zero bearing clearances shows higher bearing load peaks 
than with the nonlinear system. In the linear case, the load 
peak occurs, as expected, at the critical speed. For the 
nonlinear cases other results, not shown here, show that the 
bearing forces at bearing 2 are consistently higher than those 
at bearing 1. 

Figs. 10-a and 10-b show the Y-displacements of the rotor 
relative to the housing at bearing 1 for the start-up condi- 
tions for the linear and nonlinear cases, respectively. 

4 . 2 Periodic Response 

The harmonic balance/FFT method was used to generate 
forced nonlinear periodic responses of the generic model (Fig. 
5) in the range of 0-40,000 rpm. Figure 11 shows the funda- 
mental synchronous responses at the two critical speeds to 
occur around 1900 rad/sec an 4300 rad/sec. The figure shows 
the radial displacements of the rotor relative to the housing 
at the left bearing. The critical speed map is depicted in 
Figure 12 for the generic model in absence of bearing clear- 
ances. The gyroscopic terms are included which are shown to 
raise the forward critical speeds and lower the backward criti- 
cal speeds. The figure shows the first two critical speeds to 
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Figure 11. Critical speed of linear generic model at bearing 
1; eccentricity = 0.5 mils, no damping, no 
gyroscopic effects. 
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be exactly the same as those of the linear case obtained using 
the harmonic balance/FFT method. Numerical calculation is 
performed to investigate the existence of subharmonic 
responses . 

The trajectory of the shaft center of a 1/n order subhar- 
monic response lies in the region bounded by circles with radii 
which are determined by the harmonic and subharmonic ampli- 
tudes. The trajectory touches the bounded circles at (n-1) 
points for forward whirl and (n+1) points for backward whirl 
(See Tondl, 1973). The shaft center trajectory is studied to 
identify the subharmonic responses and to compare the amplitude 
between the harmonic response and the 1/n order subharmonic 
response with various values of the side force, gap size and 
eccentricity near the critical speeds. To allow the occurrence 
of the subharmonic response, the damping is set to a small 
value for all cases. The results are shown in Figures 13, 14, 
15 and 16. These figures show that all subharmonic responses 
occur around the second critical speed and that all are of the 
1/2 order. In Figure 13, the harmonic and subharmonic ampli- 
tudes are almost of the same magnitude for a small side force. 
In Figure 14, the harmonic response becomes dominant when 
reducing the bearing clearance size. On the contrary, the 1/2 
order subharmonic response becomes larger than the amplitude of 
the harmonic component by increasing the clearance size as can 
be seen in Figure 15. The gap size has a significant effect on 
the existence of the 1/2 order subharmonics in the HPOTP model 


Z AXIS [IN] 
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Figure 13. Trajectory for subharmonic motion of order 1/2; 
gap = 0.5 mils, eccentricity = 0.5 mils, 
side force = 10 lbs., speed = 4210 rad/sec. 
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Trajectory for subharmonic motion of order 1/2; 
gap = 0.2 mils, eccentricity = 0.5 mils, 
side force = 10 lbs., speed = 4210 rad/sec. 
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Figure 15. Trajectory for subharmonic motion of order 1/2; 
gap = 0.7 mils, eccentricity = 0.5 mils, 
side force = 10 lbs., speed = 4210 rad/sec. 
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Figure 16. Trajectory for subharmonic motion of order 1/2? 
gap = 0.5 mils, eccentricity = 0.5 mils, 
side force = 100 lbs., speed = 4210 rad/sec. 
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regardless of the side force. In Figure 16, the side force is 
raised to ten times that of Figure 13. In that case, the 
results show the amplitude of the subharmonic response becoming 
slightly larger than the harmonic response. This result 
reveals that the side force also has some influence on whether 
any subharmonic response occurs. 

More detailed investigation comparing the effects of the 
side force and gap size is studied and the results are shown in 
Figure 17. The figure shows that the side force has a signifi- 
cant effect on the existence of the 1/2 order subharmonic 
response. Around 200 lbs. side force and 0.5-0. 7 mils of gap 
size is the region which induces the large amplitude of the 1/2 
order subharmonic. Excessive side force, larger than 400 lbs. 
will cause a reduction of the subharmonic responses. With 
constant gap size and eccentricity, there exists a certain 
range of side force region in which larger subharmonic response 
would occur. This agrees with results presented by Choi and 
Noah (1987). 

The effect of eccentricity on the existence of the 1/2 
order subharmonic response with variation of the side force is 
shown in Figure 18. The figure shows that the eccentricity has 
the effect of shifting the location of the maximum subharmonic 
response to the right. The side force ranges at which the 
maximum 1/2 order subharmonic response occurs are different for 
each eccentricity. The figure also shows that the larger 
eccentricity is the lower the subharmonic response will be. 
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Gap size effect on subharmonic motion of order 1/2; 
eccentricity = 0.5 mils, speed = 4210 rad/sec. 
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This may be explained by the fact that larger response 
amplitude due to larger eccentricity induces more contact 
causing the system to become relatively weakly nonlinear 
resulting in a lesser subharmonic amplitude. It may therefore 
be necessary to reduce the gap size and side forces in any 
proposed design or maintainance criteria in order to eliminate 
dangerous subharmonic responses. 
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5. CONCLUSIONS AND RECOMMENDATIONS 


5.1 Conclusions 


The hybrid convolution approach developed in this study is 
shown to provide an efficient closed form integral formulation 
for determining the transient response of linear systems 
coupled through local nonlinearities. A typical application in 
which the present method proved quite effective is the determi- 
nation of the transient response of a generic model of the high 
pressure oxygen turbopump (HPOTP) of a space shuttle main 
engine (SSME) in presence of bearing clearances, constituting 
the local nonlinearities. Substantial savings in computation 
time were achieved as compared with direct numerical integra- 
tion techniques. 

The use of the transition matrix allows the representation 
of rotors involving skew-symmetric matrices of gyroscopic loads 
or other nonconservative systems with general velocity coeffi- 
cient matrices. A Duhamel integral would represent quite 
effectively other systems with classical modes, such as the 
housing of the HPOTP or other nonrotating, proportionally 
damped structures. The convolution formulation allows accom- 
modating with ease changes in the nonlinear or linear coupling 
parameters among the various linear subsystems involved. 

Possible improvement of the method could be achieved 
through replacement or optimization of the iteration procedure 
utilized in this study. 
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The results shown here show that the presence of bearing 
clearance can drastically alter the rotor response in the 
HPOTP. As noted by Ishida £t al^. ( 1987) for a simpler system, 
it is more difficult to go through a critical speed in the 
presence of the nonlinearity. 

An effective numerical harmonic balance algorithm for 
determining the steady-state forced vibration of a highly 
nonlinear multi degree-of-freedom rotor system has been devel- 
oped. The results show that the method is computationally 
superior to that of any direct numerical integration method. 

In addition, the complicated nonlinear steady-state periodic 
motions of multi degree-of-freedom rotor systems can readily be 
studied using the present method. 

It is shown that dangerous subharmonic resonances may 
occur for the HPOTP in the presence of bearing clearances, mass 
eccentricities, and fluid side forces. The results show that 
the ranges of the above three parameters are mutually related 
in influencing the subharmonic occurrence. In particular, the 
bearing clearance size to shaft amplitudes at the bearing, the 
local stiffness at the bearings, and the side forces, were 
found to highly influence the degree of nonlinearity in the 
rotor system and therefore determine its response and 
stability . 
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5.2 Recommendations 


It is recommended that future studies should include the 
following : 

1. Represent rotor and housing models in terms of their 
linear modal coordinates as utilized in the turbo- 
pumps of the SSME. 

2. Develop computational-Floquet or alternative analysis 
techniques for determining the stability of the 
periodic synchronous and subsynchronous response of 
rotor-housing systems. Further refine the computa- 
tional harmonic balance method for achieving more 
controlled accuracy, particularly within the reso- 
nance parameter regions. 

3. Develop method for determining the aperiodic response 
and criteria for detecting onset of probable chaotic 
motions . 

4. Perform extensive parametric studies of the steady 
state (periodic, aperiodic and chaotic) and transient 
responses of the actual SSME turbopumps to completely 
characterize their dynamic behavior for existing and 
proposed design modifications. 



57 


ACKNOWLEDGEMENT 


This work was carried out as part of a research project 
supported by NASA, Marshall Flight Center, under contract No. 

• NAS8-36293. The author and co-research assistants are grateful 
to Thomas Fox, the technical monitor, for his enthusiastic 
support and interest. 



58 


REFERENCES 


Adams, M.L., 1980, "Non-Linear Dynamics of Flexible Multi-Bear- 
ing Rotors," Journal of Sound and Vibration, Vol . 71, pp. 
129-144. 

Beatty, R.F., 1985, "Differentiating Rotor Response Due to 
Radial Rubbing," ASME Journal of Vibration, Acoustics, 
Stress, and Reliability in Design , Vol. 107, pp . 151-160. 

Bently, D., 1979, "Forced Subrotative Speed Dynamic Action of 
Rotating Machinery," ASME paper No. 74-ET-16, Petroleum 
Mechanical Engineering Conference, Dallas, Texas. 

Cipra, R.J., and Uicker, J.J., Jr., 1981, "On the Dynamic 

Simulation of Large Nonlinear Mechanical Systems, Part 1: 

An Overview of the Simulation Technique. Substructuring 
and Frequency Domain Considerations," ASME J. of Mechanical 
Design , vol. 103, pp . 849-865. 

Childs, D.W., 1978, "The Space Shuttle Main Engine High Pres- 
sure Fuel Turbopump-Rotordynamic Instability Problem," ASME 
Journal of Engineering for Power , Vol. 100, pp . 48-51. 

Childs, D.W., 1982, "Fractional-Frequency Rotor Motion due to 
Nonsymmetric Clearance Effects," ASME Journal of Energy and 
Power , Vol. 104, pp . 533-541. 

Childs, D.W. and Moyer, D.S., 1984, "Vibration Characteristics 
of the HPOTP (High Pressure Oxygen Turbopump) of the SSflE 
(Space Shuttle Main Engine), ASME Paper No. 84-GT-31, 
International Gas Turbine Conference, Amsterdam, Nether- 
lands, June. 

Choi, Y.S., and Noah, S.T., 1987, "Nonlinear Steady-State 

Response of a Rotor Support System," ASME Journal of Vibra- 
tion, Acoustics, Stress, and Reliability in Design , Vol. 
109, pp. 255-261. 

Choi, Y.S. and Noah, S.T., 1988, "Forced Periodic Vibration of 
Unsymmetric Piecewise-Linear Systems," J. Sound and Vibra- 
tion , Vol. 120, No. 3. 

Clough, R.W. and Wilson, E.L., 1979, "Dynamic Analysis of 

Large Structural Systems with Local Nonlinearities," Compu - 
ter Methods in Applied Mechanics and Engineering, Vol. 
17/18, pp. 107-129. 



59 


Davis, L.B., Wolfe, E.A. , Beatty, R.F., 1984, "Housing Flexi- 
bility Effects on Rotor Stability," MSFC Advanced High 
Pressure O 2 /H 2 Technology Conf. Proceedings, G. Marshall 
Space Flight Center, Huntsville, Alabama. 

Day, W.B., 1987, "Asymptotic Expansions in Nonlinear Rotordy- 
namics," Quarterly of Appl . Math., Vol . XLIV, No. 4, pp . 
779-792. 

Dennis, J.E., Jr., and Schnabel, R.B., 1979, Quasi-Newton 

Methods for Unconstrained Nonlinear Problems, Short Course 
on Modern Computational Techniques for Nonlinear Uncon- 
strained Optimization , Madison, Wisconsin. 

Ehrich, F.F., 1966, "Subharmonic Vibrations of Rotors in Bear- 
ing Clearance," ASME paper No. 66-MD-l, Design Engineering 
Conference and Show , Chicago, 111., May 9-12. 

Fan, U.J., and Noah, S.T., "Vibration Analysis of Rotor Systems 
Using Reduced Subsystem Models," To appear in the AIAA 
Journal of Propulsion and Power . 

Glease, J.R. and Bukley, P. , 1984, "Effect of Bearing Deadbands 
on Bearing Loads and Rotor Stability," MSFC Advanced High 
Pressure O 2 /H 2 Technology Conference Proceedings, G. 
Marshall Space Flight Center; Huntsville, Alabama, June 
27-29. 

Ishida, Y. , Ikeda, T. and Yamamoto, T. , 1987, "Transient Vibra- 
tion of a Rotating Shaft with Nonlinear Spring Characteris- 
tics During Acceleration through a Major Critical Speed," 
JSME Int . J. , Vol. 30, No. 261, pp . 458-466. 

Kubomura, K., 1985, "Transient Loads Analysis by Dynamic Con- 
densation," ASME J. Applied Mechanics, Vol. 52, pp . 559- 
564. 

Levy, A. V. , Montalvo, A., Gomez, S., and Calderon, A., 1978, 
"Topics in Global Optimization," in Lecture Notes in 
Mathematics , Note #909, pp. 18-33, Springer-Verlag, Berlin 
Heidelberg, New York. 

Meirovitch, L., 1980, Computational Methods in Structural 
Dynamics , Sijthoff and Noordhoff. 

Muszynska, A., 1984, "Partial Lateral Rotor to Stator Rubs," 
Proceedings of the I. Mech. E. 3rd. Inti. Conf. on Vibra- 
tions in Rotating Machinery, Univ. of York, 11-13 Sept., 
pp. 327-335. 



60 


Nataraj, C. , and Nelson, H.D., 1987, "Periodic Solutions in 
Rotor Dynamic Systems with Nonlinear Supports: A General 

Support," ASME Design Conference, Oct. 

Neilson, R.D. , and Barr, A.D.S., 1988, "Response of Two Elas- 
tically Supported Rigid Rotors Sharing a Common Discontin- 
uously Nonlinear Support," Proceedings of the Institution 
of Mechanical Engineers, Heriot-Watt University, 13-15 
Sept., pp. 589-598. 

Nelson, H.D. , Meacham, W.L., Fleming, D.P. and Kascak, A.F., 
1982, "Nonlinear Analysis of Rotor Bearing Systems Using 
Component Mode Synthesis," ASME Paper No. 82-GT-303. 

Noah, S.T., 1984, " Rotordynamic Analysis of the SSME Turbopumps 
Using Reduced Models," Final Report on NASA Contract NAS 8- 
34505, Texas A&M University, Sept. 

Noah, S.T., 1986, "Hybrid Methods for Rotordynamic Analysis," 
Final NASA Report, G.C. Marshall Space Flight Center, 
Alabama, under contract No. NAS8-36182, December 1986. 

Noah, S.T., Chiang, I.F. and Kim, Y.B., 1988, "Dynamic Analy- 
sis of Nonlinear Rotor/Housing Systems," MSFC Advanced High 
Pressure O 2 /H 2 Technology Conf. Proceedings, G. Marshall 
Space Flight Center, Huntsville, Alabama. 

Noah, S.T., Fan, U.J., Choi, Y.-S., and Fox, T., 1986, "Effi- 
cient Transient Analysis Methods for the Space Shuttle Main 
Engine Turbopumps," NASA Conf. on Advanced Earth-to-Orbi t 
Propulsion Tech., G. Marshall Flight Center, Huntsville, 
Alabama, May 13-15. 

Nordmann, R. , 1975, "Eigenvalues and Resonance Frequency Forms 
of Turborotors with Sleeve Bearings Crank Excitation, 
External and Internal Damping," Machine Dynamics Group, 
Technical University Darmstadt, West Germany, June 1975. 

Saito, S., 1985, "Calculation of Nonlinear Unbalance Response 
of Horizontal Jeffcott Rotors Supported by Ball Bearings 
with Radial Clearances," ASME paper No. 85-DET-33. 

Shaw, S. and Holmes, P.J., 1983, "A Periodically Forced Piece- 
wise Linear Oscillator," J . Sound Vib . , Vol . 90, pp. 
129-155. 

Tondl, A., 1965, Some Problems of Rotor Dynamics , Chapman and 
Hall, London, 1965. 

Tondl, A., 1973, "Notes on the Identification of Subharmonic 
Resonances of Rotors," Journal of Sound and Vibration , 

Vol. 31, No. 1, pp. 119-127. 



61 


Tongue, B.H. and Dowell, E.H., 1983, "Component Mode Analysis 
of Nonlinear, Nonconservative Systems," ASME Journal of 
Appl . Mech. , Vol. 50, pp. 204-209. 

Von Pragenau, G.L., 1981, "Large Step Integration for Linear 
Dynamic Systems," Conference Proc. IEEE Southeastern '81, 
reprint, April. 

Yamamoto, T.T., 1954, "On Critical Speeds of a Shaft," Memoirs 
of the Faculty of Engineering, Nagoya (Japan) University, 
Vol. 6, No. 2. 

Zalik, R.A., 1987, "The Jeffcott Equations in Nonlinear Rotor- 
dynamics," NASA/ASEE Faculty Fellowship Report, NASA 
Marshall, Huntsville, Alabama. 


62 


APPENDIX A 

DISCRETIZED FORM OF TRANSITION MATRIX FORMULATION 


The solution of the first order differential equations in 
vector form, or 

{u } = [ a] {u } + {P } (Al ) 

can be written as 

{u ( t ) } = e [a]t {u(0) } + 0 f t e [al(t " T) {PC t) } dx (A2) 

The solution of Eq . (A2) is casted in discretized form as 
{u(t. +1 ) } = e [ot]T { u(t.)} + / ^ + 1 e [a](t i+1 -t) | p(t) j dt# 

(A3) 

where T = ti + i - ti 

Take the forcing function as linear within each increment T, or 

t-t • 

{P<t)}= {P(t i )}+ ({P(t. +1 )}~ {P(t.)}) (A4) 

Let [Q ( t ) ] = e C a]T 

Substitute (A4) into (A3) to yield 

( u(t i+l ) } = tO( T )Hu(t.) } + / t . +1 e C ^ (ti + l“ t> 

t-t. 

[|P(t.) + — ^ ({P(t. +1 ) - {P(t.)})] dt 
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Then 


= [Q(T) ] {u(t. ) } + (/ t J +1 e [a] ( t i + l" t) dt ){P(t.) } 


( A5 ) 


ta]( t i + l~ t ) * [a],T- [a](t-t.), 


. ti+1 e [c](t 1+1 -t) dt . e [cOT-[a](t-ti) dt 


(/ ti+1 e [«]T-[«](t-ti) d[a]t) . U] -1 


(J tl+1 ti) d ( [o]T . [a](t-ti) })[a] _1 

ti 


= _ _ [a]T-[a] (t-ti) 


fc i+l _i 
[a] 1 

ti 


= _ ^ e [a]T-[a] (ti + i-ti) _ e [a]T-[a] (ti-ti) -j . [a] -l 


(e [0] - e [o]T ) [a ]" 1 


(e [a]T - ['I J ) [a] 1 


= (C(T) - ['I J ) [a] 


-1 


(A6) 
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Similarly , 

/ i+1 e t a Hti + x-t) - j i+1 e [a ^ T-[a] (t-t. )dt 

ti ti 1 


e [«]T ^i-l e [a](ti-t) (ti . t) d (t .. t) 
ti 


t . 


= e 


[a]T ^ i+1 e [a] (t i _t) [a](ti-t) d { [a ] ( t A -t ) }) [a ] " 2 
ti 


= e 


[a] T (e -[a]T{_ [a]T _ riv] | + [ v I v ] ) [a ] ” 2 


- [a ] -1 T + (C(T) - ['!.]) [a] -2 


(A7) 


Substituting (A6) and (A7) into (A5) yields 


|u(t. +1 )} = [C( T) ] { u ( t i ) } + ( [Q(T) ] - [ v I.])[a]" 1 ({P(t i )} + 


‘ { P(t . . )} - { P(t . )} , 

[a]" 1 y — ) - [a] _1 ({P(t i + 1 )} - { P( t ± )} ) 


( A8 ) 
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APPENDIX B 

DISCRETIZED FORM OF THE DUHAMEL INTEGRAL SOLUTION 

Consider the uncoupled set of differential equations in 
terms of the generalized coordinates {q } 


{q} + ID1 tq } + PA.] {q } = (P} 


(Bl) 


where 


[D] = PC.] = ['2Co> n .] 


['AJ = [Vj 


<u n : undamped natural frequency 


The solution of (Bl) is 


{q ( t ) } = [ 'e ^ nt cos ^dt'HqP)} + t'-rp e ?Wnt sin u> d tj» 


Uq(0) } + €co n {q(0) }) 


+ / [ v -j e C< " n(t T) sin w d (t-t) .] {p( x) }dx 

0 <3 


( B2 ) 


where = damped natural frequency. 


Let 


{h( t ) } = / ['-1 e Un(t ' T) sin w d (t-x).]{p(x)} 


‘"d 


dT 


(B3) 


P(t) is linearlized so that 


T-t 


iP(t)| - {Pj + -ji (P 1 + 1 - Pi)} 
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Substituting P(t) into (B3) and casting it in discretized form, 
one obtains 


{h(t)} = { h( t . ) } + / tl + 1 r-± e ^ {t ~ T) sin Ud (t-x)J* 

ti wd 

^ (P i + l “ p i>^ 
or 

{h( t ) } = {h(ti)} + {Ah(t)} 

For the velocity, 

|h(t)} = { h ( t i ) } + {a h( t ) } 


(B3) 


(B4) 


(B5) 


Let h(t) and Ah(t) stand for elements in vectors { h ( t ) } and 
{ Ah ( t ) } , respectively, and for an element of { P^ } , then 


Ah(t) 


_1 

od 


/ tl+1 e -«“ntt-T) sln ud(t . T) . 
ti 



X-t 


(P; 


i + 1 


P i ) ]dx 


1 -5w n (t-x) r fc i + l £u> n x • ti. \ 

= — e n j e n sin ^(t-x) • 

40(3 ti 

x- 1 . 

[Pi + <P i + 1 - Pi)]dx ( B6 ) 

Ah ( t ) = — e“ Ca,nt (INTI + INT2 ) (B7) 

a>d 
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In INTI , 

/ ^ + ^ e^ Wnt sin <D(j(t-T) dx 
t i 

= / e^ Wnt (sin cos - cos w^t sin ) dx (B 8 ) 

ti 

fc i+l a T ti +1 _ 

« sinbt / e a cos bxdx - cos bt / e (sin bx)dx (B 9 ) 

ti fc i 

where a = 5 w n and b = 

Then, from eq. (B 6 ),(B 7 ) and using equation (B 9 ) (BIO) 


INTI 




p i)] 


/ * ^e^ WnT sin uj^ (t-x) dx 
ti 


= [P i - ^ (P i+1 - p.) ] { s * nbb [e ati+1 (a cosbti+i + b sinbti+i) 


- e at i(a cosbt^ + b sinbti)] 


c ^ sb ^ [e ati+1 (a sinbti + l + b cosbti+i) 
a z +b^ 


- e afc Ma sinbti - b cosbti)]} 


In INT 2 , 

/ * + 1 Te^ nT sin (t-x) dx 
ti 


= / i + ^ ie^ nT (sin o>(j t cos 103 x - cos u),jt sin w^x) dx 

ti 


(Bll) 
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sinbt j xe ax cosbx dx - cosbt f xe ax sinbxdx 


at 


- - ^ nb -| { Te at (a cos b T + b sin b T ) ^ ~ • 

a 2 +b 2 a 2 +b^ 


n 9 

(a cos bx + 2ab sin bx ~ b cos bx ) } 


ti+1 

t i 


~ C f St> 2 { Teat ( a sin b x - b cos bx) 1 2 


a +b 


a +b 


2 2 
(a sin bx - 2ab cos bx - b cos bx)} 


b i+1 


From eq. (B6),(B7) and rearranging 


INT2 = — “ / xe 5wnX sin ud(t- T )dx 


P — P 

i+1 i ,sinbt 


e 


2,2 i+1 

a +b 


at i + l(a cos bti + 1 + b sin bti+i) 


^ 2 . -4- 2. p 2 

j~ (a cos bti + i + 2ab sin bti + i - b cos bti + i) 

a Z +b Z 


at i 

ti e ati (a cos bti + b sin bti) + ^ 

a z +b 


2 , ,2 t ,, cosbt 

(a cos bti + 2ab sin bt ± “ b cos bt i)J 2 2 

a +b 


[t i+l e 


at i + l 


(B12) 


(a sin bti+x “ b cos bti+i) 
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i ati + 1 , 2 . 


2 2 
a^ + t> 


(a sin bt i+1 - 2 ab cos bt 


2 


i +1 • b sin bt i + i> 


at 


at. 


- t^ e 1 (a sin bti - b cos bti) + 


2 2 
a +b^ 


o 

(a sin bti - 2 ab cos bti - b z sin bt^ ) ]| 


( B13 ) 


The sum INTI + INT2 is calculated from equations (Bll) and { B 1 3 ) . 
Let 


INTI + INT2 * 


sin bt cos bt 

x x — • bUnS — 

(a^+b^) 


(a 2 +b 2 ) 


. SUMC ( B14 ) 


where 


t . 


SUMS 


- [ p i ~ — | ^ P i+1 “ P i ) 1 [e at i + l (a cos bt ^+2 + b sin bti+x) 


at 


- e 1 (a cos btx + b sin bti)] + 


,p i+r p i ) 


[t i+1 e ati+1 (a cos bt i+1 + b sin bt i+1 ) 


e a *-i+l 2 o 

—7. — 5 — (a cos bt . . + 2 ab sin bt.., - b cos bt. .) 
a 2 +b 2 1+1 1+1 1+1 


ah- ati 

t^ e 1 (a cos bt^ + b sin bt^) + — ^ ^ 

a z +b z 


2 2 

(a cos bt^ + 2 ab sin bt^ - b cos bt^)] 


= e 


at i+l 


(a cos bt^ + ^ + b sin bt^ + ^) 
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I 


SUMC 



P . + 1~P . 

+ 1+ T — ~ —5 — 5 (a 2 cos bt. ,, + 2ab sin bt.^, 

a +b 11 1+1 


- b 2 cos bt i+1 ) (-e ati+ l) 


at 


- P^ e 1 (a cos bt^ + b sin bt^) + 


p i+r p i 


a 2 + b 2 


(a* cos bt^ + 2 ab sin bt i - b^ cos bt^ e at i (B 15 ) 


t. 


“ “ f ^ p i + i “ p i)^ [e ati + l (a sin bti + j - b cos bti + j) 


( P • — p . ) 

e at * (a sin bt^ - b cos bt^)] + — — — 


[t i+1 e at i+l (a sin bt i+1 - b cos bt. +1 ) 

e at i+l 2 0 

-5 — 5- (a sin bt. - 2ab cos bt. . - b sin bt . . ) 

a 2 +b 2 1+1 1+1 1+1 


at- 


at ■ 


t. e“ v ' 1 (a sin bt. - b cos bt.) + ^ - 5—5 

1 a 2 +b 2 


2 9 

(a sin bt. - 2ab cos bt. - b sin bt . ) ] 
1 1 1 


!ati+1 (a sin bt i+1 - b cos bt i+1 ) . 
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P i + 1 


' . b i (p i+r p i ) , ‘i+i (p i + i- p i» 

*■ * 1 T 1 T> J 


P . -P 
i +1 1 


1 2 

— 2 — j sin bt i + l “ 2ab cos bt i+l 

a +b 


- b 2 sin bt . , , ) (-e at i +1 ) 


i +1 


P . -P . 

- P A e ati (a sin b^ - b cos bt^ + 1+ * — - ~^~~2 

a +b 


(a 2 sin b^ - 2ab cos b^ - b 2 sin bt^) e at i (B16) 


Substitute INTI and INT2 in eq. (B15) and (B16) into eq. (B7) 
and recall that a = S n <u n and b = one arrives at 


th(t) = 


>-at .at i+1 


e ~ e ~- L ~ rx sin bt 
b(a 2 +b 2 ) 


A i +1 


, -at p ati+i 


cos bt 


b(a 2 +b 2 ) 


B i +1 


( B17 ) 


where 


A i+i = p i+i (a cos bt i+i + b sin bt i+i> 


P i +1 ~ P i 2 2 

x — 5 — (a cos bt . . + 2 ab sin bt . . - b cos bt . , ) 

T(a^+b^) 11 11 11 


P . - P . 

+ - *4 e~ aT (a 2 cos bt.+ 2 ab sin bt . - b 2 cos bt . ) 

T ( a z +b^ ) 1 1 1 


— aT 

- P^ e (a cos bt^ + b sin bt^) 


(B18) 
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B ifl ' P i + 1 (a sin bt i + l ‘ b cos bt i + l> 


P i+1 “ 2 2 

- = 5 — (a sin bt. , - 2ab cos bt. , - b sin bt . . ) 

T( a Z +b^ ) 11 11 1 

P . - P . 

+ — = — - e” aT (a^sin bt.- 2ab cos bt.- b^sin bt . ) 

T(aV) 1 


— aT 

- P^ e (a sin bt^ - b cos bt^) 


(B19) 


Take derivative of (B17), 


• e at i+l —at -at 

Ah(t) = 2 — 5 ^ " ae sin bt + b e cos bt) A^ +1 

b(a^+b ; 1 


~ a ^-i + l —at —at 

+ = — 0 (ae cos bt + b e sin bt) B. , (B20) 

b(a 2 +b^) • 11 

For t - NT, from (B3), N is an arbitrary integer, 


h (t) = f [~ e “5«n(t-x) sin wd (t-x) P( T )>] dx 

t=NT 0 ^ 


= / (*] dx 

0 


T 2T NT 

= / [•] d x + / [•] dx + ... + / [•] dx 

0 T (N-l)T 


( B21 ) 


Let fs (t) = ®I^sinbt f fc(t ) = e~ at cosb t 
b(a z +b z ) b ( a^+b z ) 


(B22) 
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Substitute eq. (B17), (B22) into (B21) yields 


h(t) = fs(t) e aT A, - fc(t) e aT B, 

t=NT 1 


+ f s ( t ) e a(2T) A 2 - f c ( t ) e a(2T) B 2 


+ fs(t)e aNT A n - f c ( t ) e aNT B n 


= f s ( t ) [e aT A a + e a2T A 2 + ... e aNT A N ] 


+ f c ( t ) [e aT B 1 + e a2T B 2 + ... + e aNT B^] (B23) 

—at — aNT 

Since e = e in (B22) for t = NT, then ea. (B23) becomes 


h(t)._ NT = - s -inbt [e -a(N-l)T A + e -a(N-2)T A + . 

b(a +b ) 1 2 


. + a n ] 


+ -CS S -bt.. [e -a(N-l)T Bi + e -a(N-2)T ^ + _ + Bj 


b(a 2 +b 2 ) 


N 

( B24 ) 


. [EA ] + 
b( a 2 +b 2 ) N 


• 1EB n ] 
b(a^+b^ ) N 


(B25) 


Now, consider t = (N+1)T, equation (B21) is written as 


h( t) t=(N+l)T 


T 2T 

= / [*] d T + / [•] dT + + 

0 


NT (N+l ) T 

/ [•] d x + / [.] dt 

(N-l)T NT 


( B26 ) 
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Equation (B23) is then written as 

h(t) t=(N + l)T “ £s(t) te3T fl l + + e3NT a m + e a(N+1)T A n+1 1 


, c /i.x r aT „ , aNT „ , a(N+l)T _ , 

+ f c ( t ) [e B 1 + + e B n + e B N+1-* ( B27 ) 

Since e at = e a ( N+ i)' 1 ’ # e q. (B26) takes the form 


h( t) t=(N+l)T 


sin bt f -aNT -a(N-l)T 

k/ 2. le A + e A + . 

b(a +b ) 


* + a n+i ] 


cos bt r _-aNT „ , „-a ( N-l )T ^ „ i 

bU^bh 1 2 *” N+1 


+ 


sin 

bt 

b(a 2 

+b 2 ) 

cos 

bt 

b ( a 

2 2 
+b 

sin 

bt 

2 

b(a z 

+b 2 ) 

cos 

bt 

2 

b(a^ 

+b 2 ) 


[e (e 


B 1 + e ~ a(N “ 2)T B 2 + . 


[e 


[e 


ea n + a n+i j 


eb n + EB N+1 ] 


• a n^ + a n+i^ 


• b n ) + B n+1^ 


(B28) 


From (B25) and B28) the following equation for t=(N+2)T can be 
written as 


h ( t ) = a + A 1 

mt, t=(N+2)T . , 2,2 . te ^N+l + A N+2 J 

b( a +b ) 
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cos bt r -aT 

+ 2 2~ e * ^N+l + ^N+2 

b( a z +b^) N 1 w ^ 


where 


ea n+i ” e ' ea n + A N+1 


EE W * e ' aT • eb n + Vl 


and EA N , EB n are defined in eq. (B24) and (B25), they are: 


. .-«(N-1)T + e -a(N-2)T + 

N 1 2 N 


EB = e - a < N “ 1 > T b + -a(N-2)T B + B 

N 1 2 N 


And A , B n are defined in eq . (B18), (B19). 


Finally, one arrives at 


h(t> ’ ■ — • ea h + 53 


. EB 


b{ a +b ) 


b( a +b ) 


N 


; t = NT, T is the time step 


where EA N , EB^ are defined in (B31) and can be calculated 
directly from previous EA N _^ and EB N _^ of time t = (N-l)T 
that is, 

EA m = e~ aT • EA m . + A m 
N N-l N 


eb n = 


-aT 


EB 


N-l 


+ B 


N 


( B29 ) 


( B30 ) 


(B31 ) 


( B32 ) 


(B33) 
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The same procedure can be followed to obtain h(t) starting from 
eg. (B20) with the result 


h(t) = 


-a sin bt + b cos bt 
b(a 2 +b 2 ) 


, a cos bt + b sin bt pR 
N b( a 2 +b 2 ) * N 


; t = NT, T is time step 


( B34 ) 


Eq. (B32) is written in vector form and is substituted into 
( B2 ) to yield 


{q( t ) } = re'^ nt cos w d tj{q(0)} + e“ €Wnt sin u> d tj • 

C {q C 0 ) } + 5w n {q ( 0 ) } ) 

+ { Si - y b v~- . EA + C9 - | - -y - . EB } t = NT (B35) 

b(a 2 +b 2 ) N b( a 2 +b 2 ) 

For {q}, derivative of(B2) for the first two terms at the right 
hand side is taken and then substitution is made of eq . (B34) 
into the resulting relationship to yield 

{q ( t ) } = p~e - ^ Wnt ( £<u n cos u» d t + w d sin o> d tj{q(0)} 

+ (- Cw n sin ui d t + cos w d t)_] . 

( {q( 0 ) } + U n {q( 0 ) } ) 


{ 


-a sin bt + b cos bt 
_ - — 

b ( a^+b z ) 


EA 


a cos bt + b sin bt 
N b ( a 2 +b 2 ) 



; t = NT 


( B36 ) 
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Equations (B35) and (B36) are written in discretized form, and 
recalling that a = S n w n , b = <*>< 3 , t = t^+i = (N+1)T, one can 
write 

{<3 ( i + l ) } = Pe at i+l cos bt i+1 ..] {q( 0 ) } 


+ P £ e ati+1 sin bt i+1 J ( { q ( 0 ) } + a{q( 0 )}) 


sin bt. , 

+ t — 2— r 1 • ea 
b(a z +b^) 


N+l 


cos bt. , 

+ — • EB H + li 


b(a 2 +b 2 ) 


(B37) 


{q( fc i + i ) 1 = p~ e atl + 1 (a cos bt. +1 + b sin bt . J _ 1 . ] {q ( 0 ) } 


i+1 


’i+1 


+ [ v e at i + l (- A a sin bt^ + 1 + cos bt v ] ( {q ( 0 ) } + 


a {q ( 0 ) } ) 


-a sin bt . ,, + b cos bt.^, 

+ t ~ ■ EA Mi , 

b(a 2 +b 2 ) N+1 


a cos bt . + b sin bt.,, 

l+l l+l 

2 2 
b(a z +b z ) 


* eb N+1^ 


( B38 ) 


where 


EA, EB are defined in (B30),(B31) and A, B in (B30),(B31) are 
defined in (B18), (B19). 
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APPENDIX C 

FOURTH ORDER RUNGE-KUTTA WITH ITERATION 

The equations ofmotion for both housing and motor can be 
cast in the first order form 


= f(t,y) 


(Cl) 


For the coefficients K 2 and K 3 of the 4th order Runge- 
Kutta Method, the coupling forces in the right hand side of 
equation (Cl) are assumed to be linear within each time incre- 
ment T = t£ + i - tf. The physical coupling forces, {Fj} at 
time t^ + T/2 are assumed then as 


i F i (t i + 1 >1 ■ -I + F i< t i + i>) 


